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For exchange-correlation functionals that depend explicitly on the Kohn-Sham orbitals, the po- 
tential V^co-(r) must be obtained as the solution of the optimized effective potential (OEP) integral 
equation. This is very demanding and has limited the use of orbital functionals. We demonstrate 
that instead the OEP can be obtained iteratively by solving the partial differential equations for the 
orbital shifts that exactify the Krieger-Li-Iafrate (KLI) approximation. Unoccupied orbitals do not 
need to be calculated. Accuracy and efficiency of the method are shown for atoms and clusters using 
the exact exchange energy. Counter-intuitive asymptotic limits of the exact OEP are presented. 
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Density functional theory (DFT) has become one of 
the most successful methods for electronic structure cal- 
culations. This success is largely due to the availabil- 
ity of a "Jacob's ladder" Q of approximations to the 
exchange-correlation (xc) energy functional E xc . High 
hopes for further improvements in functional accuracy 
rest on third-generation functionals, i.e. .functionals that 
include full or partial exact exchange 0, : One more 
exact contribution to the total energy can be taken into 
account, and long-standing DFT problems like the self- 
interaction error and the inexact high-density limit can 
be solved. Two problems so far stand in the way of a 
widespread use of the exact exchange energy in DFT cal- 
culations. Finding a correlation functional that is com- 
patible with exact exchange is the first; compatibility 
can be a severe problem for the atomization energies of 
molecules. Hyper-generalized gradient approximations 
and hybrid functionals are promising approaches to 
solve this problem on either a non-empirical or empirical 
level. But in order to use orbital-dependent functionals 
in self-consistent Kohn-Sham calculations, a second prob- 
lem - constructing the corresponding (spin-dependent) 
exchange-correlation potential V XC(T (r) - must be solved 
as well. In fact, solving the second problem might even 
help with the first, since known constraints on the poten- 
tial might be used to constrain the form of E c . But con- 
structing Vxco-(r) for a functional that depends explicitly 
on the Kohn-Sham orbitals (and thus only implicitly on 
the density) is not straightforward: The potential cannot 
be calculated directly as the functional derivative of the 
exchange-correlation energy with respect to the density, 
but must be obtained from the OEP integral equation 
0, 0, A direct solution of this integral equation is 
very demanding and so far has only been achieved for ef- 
fectively one-dimensional systems with spherical symme- 
try 0, H, 0,0- ^ or son d state systems, implementations 
within the atomic sphere approximation again reduce the 
OEP problem to spherical symmetry Q. For systems 
without spherical symmetry, the OEP has recently been 



constructed by directly evaluating the response function 
lill LUJ, list LUI ■ This makes molecular calculations possi- 
ble, but requires not only occupied but also unoccupied 
orbitals and may lead to instabilities in the asymptotic 
regions of the OEP, as discussed in [13. Hjl Hq , so al- 
ternatives continue to be developed |l7|. On the other 
hand, the KLI and localized Hartree-Fock [lil lisj ap- 
proximations are easier to implement than the full OEP, 
yet accurate in many cases [(J |3, |la |l8|, |l9|, |20( . 

We demonstrate how V xcrT (r) = SE XC /Sn a (r) can be 
calculated for any orbital-dependent E xc without explic- 
itly solving an integral equation. The exact OEP is ob- 
tained iteratively by self-consistently solving a system 
of partial differential equations that is coupled to the 
Kohn-Sham equations. The iteration converges quickly. 
We use an essentially asymptotically-correct approxima- 
tion as the starting point, and the correct asymptotic 
behavior is preserved during the iteration. The unoccu- 
pied Kohn-Sham orbitals need not be calculated. The 
orbital shifts we evaluate were defined in Refs. 0, ||, 
but have not been calculated before. However, orbital 
shifts have been used successfully in linear response the- 
ory [2ll , and the advantages of avoiding explicit eval- 
uation of the response function are discussed in Ref. • 
We present self-consistent exact-exchange OEP calcula- 
tions for atoms and three-dimensional sodium clusters to 
demonstrate the accuracy and efficiency of the method. 

Excellent reviews of the OEP method have been pub- 
lished 0, |2(J , so we restrict ourselves to the central equa- 
tions. Following 0,0, the OEP integral equation is writ- 
ten in the form 

$>t,(r) W «,(r) + c.c. = 0. (1) 

»=i 

Here, c.c. denotes the complex conjugate of the previous 
term, and (^ia-(r) denotes Kohn-Sham orbitals, i.e., the 
solutions of 

(h K Sa - Sic) Pia(r) = 0, (2) 
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where /ikSct — — (^ 2 /2m)V 2 + V a (r) is the Kohn-Sham 
Hamiltonian, V^r) = V ex t(r) + Va(r) + V xca (r) is the 
Kohn-Sham potential with V oxt (r) the external potential 
and Vfj(r) = J d 3 r'e 2 n(r')/\r — r'| the Hartree potential. 
The unnormalized orbital shifts i/j* a (r) are defined by 



(3) 



y J ^) [WrQ - "xci. (rQ] y Jg (r') dV . 

3'=1 



where 



,(r) = 



1 <5£ xc [{^ T }] 



(4) 



(Ref. uses K(J (r) = -CWM,.) The -^(r) 
are the first-order perturbation theory corrections to the 
(p l(T (r) for the perturbing potential (r) - Vxco-(r). 
Thus, Eq. asserts that the density shift vanishes to 
first order, and perturbation theory shows that ip* a (r) 
must fulfill the partial differential equation 



-[V XC a(r) - u xclrT (r) - (V x< 



U xcla )}if* a {Y). (5) 



In Eq. |JSJ|, the orbital averages are V XC i a = 
J ^* a (r)V XC!7 (r)ip ia (r)d 3 r and u xcla = 
J tp* cr (r)u xcilT (r)if ilJ (r)d 3 r. Multiplying Eq. |TJ by 
V a (r) and inserting Eq. JSJ), solved for V„(r), into the 
resulting expression, an equation is obtained that can 
be solved for V xca (r), and a few further manipulations 
la, Si then lead to the form 



V xca (r) 



i N " r 



2n CT (r) frf 



XCMT "-XCiCT ) 



(6) 



The KLI potential is obtained by setting ip* a (Y) = OVi 
in the above equation, a choice that can be motivated in 
the sense of a mean- field approximation 0, . 

Combining Eqs. © and © to calculate V XC(7 (r) in a 
self-consistent iteration is the basic idea of our method. 
Calculating the ipia( r ) 1S thus the first crucial step. The 
first problem to be addressed is which boundary condi- 
tions are appropriate. Since the VCr( r ) are orbital shifts, 
the boundary condition to be employed in the solution of 
Eq. JSJ must be the same as the one used in solving Eq. 
(J2J), i.e, for our finite systems, VCr( r ) — > for r — > oo |8|. 
But at first sight, calculating ip* a {r) still seems a hopeless 
task: Evaluating Eq. is highly impractical and Eq. (JSJ 
is singular. It does not have a unique solution because to 
any one solution ipi a (r), the Kohn-Sham orbital <Pi a {r) 
multiplied by an arbitrary constant can be added to find 
another solution. 



That Eq. J5J) can nevertheless be solved efficiently is 
based on two facts. First, it is clear from Eq. that 
the particular solution of Eq. JSJ that is relevant for 
the OEP must be orthogonal to the Kohn-Sham orbital 
ipi a (r). Second, the right hand side of Eq. (JSJ is or- 
thogonal to ifi a (r). The orbital shift therefore can be 
calculated elegantly by the conjugate gradient method 
23]. In short, if the operator on the left hand side of Eq. 
|J5J is abbreviated by Aj and the right hand side by bi, 
then the conjugate gradient method constructs the solu- 
tion of Aiip* = b{ by calculating a sequence of additive 
corrections to a starting guess ip* s . The first correction 
in the sequence is, up to a constant factor, given by the 
expression 6, — Aiip* s . Since bi is orthogonal to ifi, and 
since tpi spans the nullspace of Ai, the whole correction 
is orthogonal to ipi. Since the following steps in the se- 
quence proceed analogously [23|, the desired orthogonal 
solution is obtained if ip* s is chosen orthogonal. In prac- 
tice, numerical inaccuracies might lead to small contri- 
butions of ifi to the calculated ipi , but they can easily be 
removed by orthogonalizing the final to if i (or, in the 
case of degeneracies, to all Kohn-Sham orbitals having 
the eigenvalue £,). Real space and basis-set implementa- 
tions should be equally unproblematic since V x {r — > oo) 
is essentially correct in the starting approximation. 

One way of constructing the OEP would be to insert 
the orbital shifts into Eq. ©. But this is numerically 
cumbersome for finite systems: Evaluating the last term 
on the right hand side of Eq. © requires dividing deriva- 
tives of exponentially decaying functions by the expo- 
nentially decaying density. For large distances where all 
orbitals have decayed to the extent that their numeri- 
cal representation comes close to zero, inaccuracies can 
be introduced into the iteration. We found that, with 
appropriate cut-offs, the iteration converged quickly for 
the Be atom. However, for non-spherical systems this 
procedure becomes cumbersome. 

An alternative way of using the orbital shifts to con- 
struct the OEP is based on the fact that once the ip* a (r) 
have been calculated for a given V xclT (r), one obtains from 
them an estimate for how close the given V xca (r) is to the 
true OEP in the region of nonvanishing density by defin- 
ing the function 



S*(t) 



E 



(7) 



For the exact OEP, Eq. Q reduces to Eq. QJ, i.e., 
S a {r) = Vr. But if fi a (r) and tp^ir) have been cal- 
culated via Eqs. an d © f° r an approximate V xca (r), 
then S a (r) does not vanish. As an example, the left half 
of Fig. ^ shows the OEP and KLI potentials for the Ne 
atom, and the full line in the right half of Fig. Q] shows 
the S(r) obtained for the KLI potential. (When there 
is no spin polarization, we drop the index a.) Although 
S(r) is not an exact representation of the difference be- 
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FIG. 1: Left: I4(r) for the Ne atom versus radial coordinate r in bohr (ao). Dashed: KLI. Full: OEP. Right: S(r) at different 
stages of the iteration. Full: S(r) for KLI potential. Long dashed, short dashed and dotted: S(r) after first, second and third 
iteration. Note that the magnitude of S(r) is reduced by each iteration and goes to zero. 



tween the two potentials, as can be seen by comparing 
S(r) to Fig. 2 of Ref. 0, it nevertheless indicates where 
and qualitatively how the approximation deviates from 
the OEP: S(r) vanishes in the asymptotic region where 
KLI is correct, but features the greater intershell bump 
and lower core value of the OEP. 

This observation motivates the following simple itera- 
tive construction of the OEP: Using an approximation to 
the exact OEP, calculate <Put(t) and V'io-OO v ^ a ^qs. @ 
and J5J and S^(r) via Eq. Q. An improved approxima- 
tion to the OEP is then constructed as 

V™(r) = V°£(r)+cS ( ,(r), (8) 

where c is a positive constant having the dimension of 
energy over density. The idea of this iteration is to re- 
move the error - for which S a (r) is a measure from 
the current approximation to obtain a better one. The 
numerical value of c is a system-dependent parameter of 
the iteration. Larger values lead to faster convergence, 
but too large a value will ultimately lead to divergence, 
c is thus similar to the parameter introduced in the over- 
relaxation method |23j or the mixing parameter in self- 
consistent calculations. In the calculations discussed be- 
low we used values of about 1 for the atoms and of about 
30 for the clusters on a trial-and-error basis. A ratio- 
nale for these values will be presented in future work. 
In a simple but important second step, it must be en- 
sured that the improved approximation satisfies the re- 
lation V^f.N^a — u XC N^ a which is obeyed by the exact 
OEP |fj, |S| • This can easily be achieved by adding the 
constant u xcN<ja - J Lp* N<ja (r)V x n °J'(r)if Niya (r)d 3 r to the 
V x ™(r) computed from Eq. JHJ. With this new approx- 
imation for V xca (r), Eqs. © and JSJ can be solved anew 
to start another cycle until convergence is achieved. 

In our calculations the most efficient way of doing one 
iteration of Eq. @ is to do several cycles of iterating 



{tp* a (r)} and V xca (r) while keeping {<p la (r)} : {u xci(7 (r)}, 
/iks an d Si in Eq. (JjJJ fixed. Only then are the Kohn- 
Sham equations solved anew. The reason for this is that 
in our real-space finite difference calculations, solving Eq. 
© requires only roughly as much computational effort 
as calculating the Hartree potential, whereas solving the 
Kohn-Sham eigenvalue equations is more demanding. 

We first employed this scheme to calculate OEP to- 
tal energies and eigenvalues of spherical atoms in the 
exchange-only approximation. The corresponding equa- 
tions were solved on a logarithmic grid. The right half of 
Fig. G] shows that each iteration reduces S, as it should. 
The first line of Table [Q shows the KLI energies, the 
starting point of our iteration. The second line gives 
the same values for the OEP, and the first column indi- 
cates how many iterations were necessary to converge to 
0.0001 hartree accuracy. Our simple iteration converges 
within a few steps to very high accuracy. We have care- 
fully checked that these numbers correspond to the ones 
obtained by the much more involved direct solution of 
the OEP integral equation, by re-running the program 
used in Refs. 0,0] on fine grids. 



TABLE I: Total energy and Kohn-Sham eigenvalues in hartree 
for the Ne (top) and Ar (bottom) atom. First line (iter=0): 
KLI approximation. Second line: OEP. The first column gives 
the iteration number at which all eigenvalues are converged 
to 0.0001 hartree accuracy. Total energies are already correct 
after the first iteration. 
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-128.5448 


-30.8021 


-1.7073 


-0.8494 
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-128.5454 


-30.8200 


-1.7181 


-0.8507 









-526.8105 


-114.4279 


-11.1820 


-8.7911 


-1.0942 


-0.5893 
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-526.8122 


-114.4522 


-11.1532 


-8.7338 


-1.0993 


-0.5908 
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TABLE II: KLI and OEP total energies and eigenvalues from 
three-dimensional cluster calculations, listed as in Table [I] 
Top: Na4, Bottom: Nag. See text for details. 



iter 
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£l 


£2 


£3 


e 4 





-0.7528 


-0.1776 


-0.1394 






5 


-0.7531 


-0.1762 


-0.1401 







-1.5282 -0.1963 -0.1591 -0.1458 -0.1458 
5 -1.5285 -0.1954 -0.1592 -0.1459 -0.1459 



and broad range of applicability of our approach. It thus 
opens a path for employing orbital functionals, in partic- 
ular the exact exchange energy, in self-consistent Kohn- 
Sham calculations. 
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Our scheme can easily be used to calculate the OEP for 
three-dimensional systems without particular symmetry. 
To demonstrate this, we have calculated the OEP for the 
sodium clusters Na4 and Nag, using the pseudopoten- 
tial and geometries from Ref. |24|. A three-dimensional 
Cartesian grid with 129 points in each direction and of 
spacing 0.5ao was employed to obtain an accuracy of 
0.0001 hartree. Eq. @ was solved by damped gradi- 
ent iteration with multigrid relaxation, and Eq. J5J) by 
our own conjugate gradient routine. In Table 1TT1 we list 
total energies and eigenvalues obtained with the exact ex- 
change functional, iterating from KLI to OEP. The first 
thing to note is that, even for the clusters, KLI is a rather 
good approximation to the OEP. But it should also be 
noted that the relative error in the KLI energies is nearly 
two orders of magnitude higher for the clusters than for 
the atoms. Nevertheless, just five steps of our iteration 
refine KLI to OEP energies. The relative error in the 
Levy-Perdew virial relation [25| falls to ~ 10~ 5 during 
the first five iterations, and to ~ 10~ 7 upon further it- 
eration, demonstrating the accuracy and stability of our 
method. 

Finally, we investigated the asymptotic behavior of 
the Kohn-Sham potential. Contrary to common belief, 
the exchange potential does not fall off asymptotically 
like —e 2 /r everywhere for finite systems, but approaches 
—e/r + C, where C is a nonvanishing constant, on nodal 
surfaces of the energetically-highest occupied orbital. For 
Na 4 we find C = 0.0307 for OEP and C = 0.0298 hartree 
for KLI. This unexpected behavior of V x (r) can be under- 
stood on the basis of Eq. It is not reproduced by any 
of the common methods such as the local density approx- 
imation, and also not by the OEP construction presented 
m Ref. hj. Our work thus demonstrates the existence 
of non- vanishing asymptotic constants in the exact OEP, 
and confirms the arguments recently put forward in Ref. 

In conclusion, our simple iterative construction of the 
OEP does not require solving an integral equation or cal- 
culating unoccupied Kohn-Sham orbitals. We discussed 
the surprising asymptotic behavior of V x (r). All-electron 
calculations for atoms and three-dimensional pseudopo- 
tential calculations for clusters demonstrate the accuracy 
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